Logarithmic Scaling of Catastrophic Collapse 
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We study the catastrophic stationary self-focusing (collapse) of laser beam in nonlinear Kerr 
media. The scaling of a self-similar solutions near collapse distance z = z c obeys a (z c — z) 1 ^ 2 scaling 
law with the well-known, leading order modification of log-log type. We show that the validity of 
the log-log modification requires double-exponentially large amplitudes of the solution ~ 10 lol °°, 
which is unrealistic to achieve in either physical experiments or numerical simulations. We derive a 
new equation for the adiabatically slow parameter which determines the system self-focusing across 
a large range of solution amplitudes. Based on this equation we develop a perturbation theory 
for scaling modifications beyond the leading log-log. We show that the new scaling agrees with 
numerical simulations beginning with amplitudes around only three times above of the initial pulse. 

PACS numbers: 42.65.Jx, 42.65.-k, 52.38.Hb 



The catastrophic collapse (self-focusing) of a high 
power laser beam has been routinely observed in non- 
linear Kerr media since the advent of lasers [3-0] ■ The 
propagation of a laser beam through the Kerr media is 
described by the nonlinear Schrodingcr equation (NLSE) 
in dimensionless form, 



(i) 



where the beam is directed along z-axis, r = (2, y) are 
the transverse coordinates, ^(r, z) is the envelope of the 

electric field, and V e The NLSE © also 

describes the dynamics of attractive Bose-Einstein con- 
densate (BEC) [5J (z is replaced by the time variable in 
that case). In addition, the NLSE emerges in numer- 
ous optical, hydrodynamic, and plasma problems, and 
describes the propagation of nonlinear waves in general 
nonlinear systems with cubic nonlinearity. 

If only one transverse coordinate is taken into account, 
then the NLSE is integrable by the inverse scattering 
transform Q leading to global existence for all solutions. 
A solution of the NLSE which depends on both transverse 
coordinates (x,y) can develop a singularity ("blow up") 
such that the amplitude of the solution reaches infinity 
in a finite distance z c . Since the blow up is accompanied 
by dramatic contraction of the spatial extent of function 
■0, it is called "wave collapse" or simply "collapse" 0,0- 
Near the singularity z — z c , the NLSE looses applicabil- 
ity, and either dissipative or non-dissipative effects must 
be taken into account. Such effects can include the opti- 
cal damage and formation of plasma in the Kerr media, 
inelastic scattering in the BEC, or plasma density deple- 
tion in high temperature laser- plasma interactions jjjllioj. 

Equation fl} can be rewritten in the Hamiltonian 
form iip t = J^r with the Hamiltonian H = j (|V^| 2 — 

^|'0| 4 )<ir. Another conserved quantity, N = J \ip\ 2 dr, has 
the meaning of the optical power (or the number of parti- 
cles in the BEC). The sufficient condition for the collapse 
is H < 0, while the necessary condition is N > N c , where 



N c is the critical power as defined below. 

While the large power N 3> N c typically produces mul- 
tiple collapses (multiple filamentation of the laser beam 
[ill ]) with strong turbulence behavior [lj], EH, the dy- 
namics of each collapsing filament is universal and can 
be considered independently. Each collapsing filament 
carries the power N only moderately above N c . We con- 
sider a single collapsing filament centered at r = 0. For 
z — >• z c the collapsing solution of the NLSE quickly ap- 
proaches the cylindrically symmetric solution, which is 
convenient to represent through the following change of 
variables H: 



if>(r,z) = t V(p,t) 



ir+iLLzp 2 /4 



Here, L(z) is the z-dependent width of solution and 

dz' 



V 



L\z>) 



(2) 



(3) 



are blow up variables such that r — > 00 as z — > z c . Trans- 
formation ([2]) was inspired by the discovery of the addi- 
tional conformal symmetry of the NLSE which is called 
the "lens transform" [14H16j |. 

It follows from Q, © and © that V(p,r) satisfies 



id T V + V 2 „V -V+ \V\ 2 V + -o 2 V = 0, 



where 



-L 3 L, 



and V 2 p = d 2 p + p^dp. 



(4) 



(5) 



As z — > z c , j3 — > adiabatically slowly and V(p) ap- 
proaches the ground state soliton R(p) [l6| . The ground 
state soliton is the radially symmetric, z-independent so- 
lution of the NLSE, -R+V 2 p R+R 3 = 0. It is positive def- 
inite, i.e., R > 0, with asymptotic R{p) = e~ p \App~ 1 / 2 + 
0(p- 3 / 2 )], p -> 00, A R ee 3.518062... 0. Also 
R defines the critical power N c = 2ir j R 2 pdp = 



2 



11.7008965 . . . The limiting behavior in V — > R as z — >• z c 
implies that the d T V term in ((4]) is a small correction 
compare to the other terms. 

Refs. 17i] and [HI found that the leading order depen- 
dence of L{z) has the following square-root-log-log form 



L 



2tt 



In | ln(z c 



1/2 



(6) 



(Ref. II 711 has a "slip of pen" in a final expression, see 
e.g. [19j for a discussion.) The validity of the scal- 
ing © at z — > z c was rigorously proven in Ref. (20| . 
However, numerous attempts to verify the modification 
of L oc (z c — z) 1 / 2 scaling have failed to give convinc- 
ing evidence of the log-log dependence (see e.g. pHI^.) 
Note that without logarithmic modification, the scaling 
(z c — z) 1 / 2 implies /3 = const, N = oo, and infinitely fast 
rotation of the phase for r — > oo with j3 ^ 0. Thus, the 
logarithmic modification is necessary and is responsible 
for the adiabatically slow approach of j3 to 0. 

A qualitatively similar problem of logarithmic modi- 



fication of square-root scaling also occurs in the Keller- 
Segel equation, which describes either the collapse of self- 
gravitating Brownian particles or the chemotactic aggre- 
gation of micro-organisms It was shown in [271 ] 
that the leading logarithmic modification in Keller-Scgcl 
equation is valid only for very large amplitudes (> 
10 10000 ). Also in [27|, the perturbation theory was de- 
veloped beyond the leading order logarithmic correction, 
which is accurate starting from moderate amplitudes 
(> 3) of collapsing solution. 

Following qualitatively some ideas of [27} , in this paper 
we develop the perturbation theory about the self-similar 
solution of (j4|) with V ~ R and show that the scaling (J6j) 
dominates only for very large amplitudes 



\ip\ > 10 



10 



100 



(7) 



Instead of pursuing this unrealistic limit, we suggest a 
new expression (derived below) as a practical choice for 
the experimental and theoretical study of self-focusing: 



r„ / mi/2 A a „, o „, , a 4(-l - 41n3 + 41nlnv4) 

L = [2tt(z c - z)] 1 ' 2 In A - 41n3 + 41nln/l + -i — — '- 

\ in J\ 

-28 - 80 In 3 - 32 (In 3) 2 - tt 2 ci + 80 In In A + 64(ln 3) In In A - 32[ln In Af 



{\nAf 



-1/2 



[2tt(z c - z)\ 



1/2 e 



L(z ) 



(8) 



M = 44.773 . . . , fa = P(z ), ci = 4.793 . . . , c 2 = 52.37 . 



a 



ev* / 2/3 2 8/3q /2 2(il (20 + tt 2 ci) 12^ /2 (20tt 3 + ir 5 a) 2/3% (840tt 3 + 42tt 5 ci + tt 7 c 2 ) ~ 



M 



r 



This expression depends on an additional parameter, zq. 
In practice, we choose zq as the smallest distance at which 
the collapsing solution has approximately reached the 
self-similar form. In numerical examples below, it implies 
that z = zq corresponds to distances where the amplitude 
of the initial Gaussian pulse has increased ^3—15 fold. 
Note that zq can be chosen arbitrarily larger than this 
minimum distance as long as it is less z c . 

To illustrate the poor agreement with the log-log law 
at moderate amplitudes, Figure [T] shows the dynamics of 
L(z) obtained from numerical simulations. The simula- 
tions were started with different initial conditions in the 
form of Gaussian beams ip(r, 0) = pe~ r . Figure Q] shows 
that L(z) neither agrees with the log-log law (J6j) nor it is 
universal. In contrast, the dependence of /3 T on /3 appears 
to be universal as demonstrated in Figure [2l The curves 
corresponding to different initial conditions converge to 
a single /3 T (/3) curve after the initial transient distance. 



These points of convergence define zq in our simulations, 
with numerical values given in the caption of Figure [1] 

To determine j3 T (/3) analytically, we consider the 
ground state soliton solution Vo(f3,p) of ^ given by 



V 2 F 



V 



V 3 + ^p 2 \n = i). 



(9) 



Vo(p,p) has an oscillating tail, 
os [^!p 2 -/?- 1/2 lnp + J +0(p-% 



The function 

V {f3,p)=cp- 1 

with c, 0o = const and p 3> 2//3 1 / 2 . Here, by ground 
state soliton Vq, we mean the real function such that 
it minimizes |c| in the tail. It implies that Vq has only 
small amplitude oscillations with |c|<§;lfor0</3^;l. 

The full solution V(/3,p) of (|4]) is well approximated 
by Vo(/3,p) for p < 1. However, the small but nonzero 
value of d T Vo = ftr^p- provides an imaginary contribu- 
tion to V. To account for the imaginary contribution at 
the leading order, we allow Vq to be complex (replacing 
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FIG. 1: (Color online) Dependence of L on z c — z obtained 
from numerical simulations of the NLSE JTJ) (solid lines) and 
from Eq. © (dashed lines) for different initial conditions. 
Each pair of closely spaced solid and dashed lines corresponds 
to the same initial condition and has relative difference ~ 2%. 
The curves are labeled by the power N for each Gaussian ini- 
tial condition xj){r, 0) = pe~ r . Values of zo in (HI are chosen 
at points where solid curves of Figure converge to the sin- 
gle universal curve and correspond to f3 = 0.20, 0.17, 0.14, 
and 0.11 for N/N c = 1.37, 1.21, 1.13, and 1.08, respectively. 
The dotted line shows L from ([6|. The inset shows L(z) for 
N/N c — 1.13 starting from the beginning of simulation, z — 0. 

it by Vq), similar to the approach of [H |lf|. We for- 
mally add an exponentially small term ii/(8)Vo to (JSJ) as 
follows, V 2 Vb - % + \V \ 2 Vo + f p 2 V ~ iv(fl)V = 0. The 
yet unknown v(B) accounts for the loss of power of Vq 
by emission into the tail. One can reinterpret the re- 
sulting equation as a linear Schrodinger equation with a 



self-consistent potential U = — \Vq\ 2 — jp 2 and a complex 
eigenvalue E = — 1 — w(B). (This type of nonself-adjoint 
boundary value problems was introduced by Gamov in 
1928 in the theory of a-decay [H)].) Assuming |3« 1, 
we identify two turning points, p a ~ 1 and pb ~ 2//3 1 / 2 , 
at which Re{E) + U = 0. Using the WKB (Wentzel- 
Kramers-Brillouin) approximation we consider the tun- 
nelling from the collapsing region p < 1 through the clas- 
sically forbidden region p a < P < Pb, and obtain, similar 
to " 



12| that 



V = e 2 £ 1/2 exp 



itL_ p 2 -0-1/2 lap- ifo 



2l/2 A 

x /4 [p- 1 + 0(p-% 4> = const, p-> Pb , (10) 

where Ar results from the matching of the asymptotic 
of R with the WKB solution. We also note that the 
tail (|10|) is derived in the adiabatic approximation which 
is valid for large but finite values of radius, 2//3 1 / 2 <C 
r/L <C A (2//3 1 / 2 ), where A{z) ^> 1 is a slowly changing 



factor in comparison with L{z). Even though for r/L > 
A(2//3 1/2 ) the solution is not self-similar fig, El Eil , its 
large-radius asymptotic has no influence of L(z) and is 
not considered here. 

We define the power (the number of particles) iV& in 
the collapsing region p < p^ as 



N h 



J |^| 2 dr = 2vr J \V\ 2 pdp. 



(11) 



r<p b L 



p<Pb 



and a flux P beyond the second turning point pi, = 2/ y/]3 
as P = 2Trp[iVV* + c.c] \ P=Pb , where c.c. stands for 
complex conjugate terms. From conservation of N , the 
flux P determines the change of as 



dNb 
dr 



-2np [iVV* + c.c] , p > p b , 



(12) 



where we approximated P at p = pb through its value at 
p S> Pb taking advantage of almost constant flux to the 
right of the second turning point. Using the adiabatic 
assumption that = Pr^ja-, and approximating V in 
fll~2l) by (OH we obtain that 



-A«Al 



dN b 
dB 



01/2 



(13) 



Recalling the definition of v(j3), one can also find v(/3) ~ 

(2nA 2 R /N b )e~^ from Q3J). 

The next step is to find ^ in (fI5]l. We based 
our derivation on a crucial observation that the abso- 
lute value |V(/3, p)| of the numerical solution of dU co- 
incides with Vo(P,p) for < p < pb, as shown in 
Figure [3l In contrast, the approximation Vq(0, p) ~ 
R(ft) + dV(fB, p)/d/3\/3=o used previously (see e.g. [16j) is 
limited to p -C pb because the amplitude c of the tail of 
Vq has the essential complex singularity c oc e'^^ 2 ^ 1 ^ 
for /3 — > 0. Approximating through replacing V in 
(fTTjl by Vo(j3,p) we obtain the following series 



dN b 
dB 



= M [l + ci/3 + c 2 /3 2 + c 3 /3 3 + c 4 /? 4 + c b B b ] , (14) 



where M = (2ir)- 1 dN b /df3 = (1/4) J °° p 3 R 2 (p)dp = 
0.55285897 ... and coefficients c x = 4.793, c 2 = 
52.37, c 3 = 296.99, c 4 = -4660.87, c 5 = 10540.4 are es- 
timated from the numerical solution of ([9]). Figure [2] 
shows that (fl3|) and (fl4|) well approximate the full nu- 
merical solution for 3 < 0.18. Indeed, B T {8) from (fl3|) 
with dNb/dB, obtained either numerically via Va{6, p) or 
using Eq. p^|) . arc indistinguishable on the plot (they are 
both shown by the dash-dotted line). Also a dashed line 
corresponds to Eq. (|14p with only c\ and c 2 taken into 
account. For comparison, dotted line shows the standard 
approximation for 8 T (0) , which corresponds to (fT4|) 
with the expression in square brackets replaced by 1. As 
we see, the standard approximation fails all way down 




FIG. 2: (Color online) Solid lines show /3 T (/3) from numerical 
simulations of the NLSE {TJ with the same initial conditions 
as in Figure [1] The curves are labeled by the values of N/N c . 
See the text for the description of dashed, dashed-dotted and 
dotted lines. The inset shows f3 as a function of z c — z for the 
same N/N c . 



to < 0.07. Here, numerical values of T (0) shown for 
> 0.07 are obtained from simulations of the growth of 
collapse up to max|-0| ~ 10 15 , and our simulations do 
not resolve the range < 0.07. 

Eqs. (|3l).(|5|). (|13|) and (fl4|) form a closed system from 
which 0{z) and L(z) can be determined. The asymptotic 
solution of this system, in the limit z — > z c , r — > oo, 
— » and L — > 0, gives our main result, Eq. ©, where 
M = 2A\/M. Deriving ©, in Eq. ((Til) we used leading 
order terms with C\ and C2 only. (The corresponding 
T {0) is shown by a dashed line in Figure [2j) Thus c\ and 
C2 in (|14l) are sufficient to produce very good agreement 
with the simulations shown in Figure [T] 

Finally, we comment on how we determine L and 
from numerical simulations. At each z with high pre- 
cision we use the following two-step procedure. First, 
we determine L(z) from the numerical solution ip(r, z) 

-1/2 



as L 



M 1 



,|V>I„ 



r=0 



an expression derived 



from the Taylor series expansion of Vo(0, p) for p <C 1 in 
©. Second, we determine /3(z) from the nonlinear con- 
dition |^>(0, z)\ — j]-zjVo{(3,0) using the pre-computed 
values of Vb(/3, 0) from the solution of ([9]). We found 
that this procedure gives much better accuracy in deter- 
mining L and than the alternative procedures reviewed 
e.g. in Ref. [16] . 

In conclusion, we found that the collapsing solu- 
tion is described by the self-similar solution \ip(r, z)\ = 
TR^o(>(^Tr) ^ < r/L(z) < 2//3 1 /2 (z) with 
L(z) given by §8§ and — —L 3 L ZZ , where Vo(f3,p) is 
the ground state soliton solution of ©. For r/L(z) ^S> 



FIG. 3: (Color online) Asymptotic p > 1 for Vo (solid line), 
full numerical solution \V\ (dotted line) and R (dashed line) 
for P = 0.073. It is seen that Vo and \V\ almost coincide for 
P < Pb- 



2/(3 1 / 2 {z) the collapsing solution has the tail Vq from 
dTUI) - We found that the dependence L(z) in © is in 
very good agreement with the direct numerical simula- 
tions, as shown in Figure [1] starting from quite moder- 
ate increase (<~ 3 times) of the amplitude of the initial 
Gaussian beam. By the direct substitution of the values 
L(z ) and (3(zo) into (JHJ) , with z defined in Figure[lJ one 
can see that expression ([5]) matches the classical result 
© with accuracy ~ 10 — 20% only for the unrealistically 
large amplitudes given by (JTJ). It suggests that the classi- 
cal result ([6]), while being asymptotically correct, should 
be replaced by much more accurate new formula ((8]) for 
any currently foreseeable physical systems. 
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